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PREFACE 


The  work  described  in  this  report  was  performed  during  the  period 
May  to  December  1995  by  Dr.  Richard  D.  Rechtien,  Consulting 
Geophysicist  and  Professor  Emeritus,  University  of  Missouri  at  Rolla.  The 
work  was  performed  under  Contract  No.  DACW39-95-M-2999  and  funded 
by  the  Earthquake  Engineering  (EQEN)  Research  Program  project, 
“Geophysical  Methods  for  Site  Characterization  and  Measurement  of 
Seismic  Properties,**  Work  Unit  No.  33014. 

This  report  was  prepared  by  Dr.  Rechtien.  Principal  Investigators  for 
the  EQEN  Research  Project  and  technical  monitors  for  this  work  were  Dr. 
Dwain  K.  Butler  and  Mr.  Robert  F.  Ballard,  Jr.,  Earthquake  Engineering 
and  Geosciences  Division  (EEGD),  Geotechnical  Laboratory  (GL),  U.S. 
Army  Engineer  Waterways  Experiment  Station  (WES).  This  work  was 
performed  under  the  general  supervision  of  Dr.  Arley  G.  Franklin,  Chief, 
EEGD,  and  Dr.  William  F.  Marcuson,  Director,  GL.  Dr.  Mary  Ellen  Hynes 
was  EQEN  Project  Manager. 

At  the  time  of  publication  of  this  report.  Director  of  WES  was  Dr. 
Robert  W.  Whalin.  Commander  was  COL  Bruce  K.  Howard,  EN. 


The  contents  of  this  report  are  not  be  used  for  advertising,  publication, 
or  promotional  purposes.  Citation  of  trade  names  does  not  constitute  an 
official  endorsement  or  approval  of  the  use  of  such  commercial  products. 


Conversion  Factors, 
Non-SI  to  SI  Units  of 
Measurement 


Non-SI  units  of  measurement  used  in  this  report  can  be  converted  to  SI  units  as  follows: 


1  Introduction 


Background  and  Scope  of  Work 


Liquefaction  is  a  term  used  to  describe  a  process  involving  the 
complete  loss  of  shear  strength  of  loose-  to  medium-dense  sand  deposits  (or 
other  non-cohesive  soils  )  below  the  water  table  during  the  passage  of 
large-amplitude  earthquake  waves.  The  existence  of  such  processes  in 
nature  is  well  evidenced  by  surface  observations  at  many  earthquake  sites 
throughout  the  world  (see,  for  example,  Kawasumi  1968,  Seed  et  al.  1990, 
and  CAEE  1995).  A  primary  objective  of  liquefaction  research  is  the 
development  of  indicators,  derived  from  in  situ  geotechnical 
measurements,  that  provide  accurate  assessment  of  liquefaction  potential  of 
a  cohensionless  deposit  in  a  given  earthquake  environment. 

Seismic  waves  of  sufficient  magnitude  to  cause  liquefaction  cannot 
(generally  speaking)  be  practically  produced  for  in  situ  liquefaction 
studies.  Consequently,  investigations  toward  identification  of  relevant 
source/soil  characteristics  have  essentially  been  confined  to  laboratory 
experimentation  on  both  field  and  synthetic  soil  samples.  While  these 
studies  have  provided  a  good  understanding  of  the  mechanisms  of  the 
stress/soil  dynamic  that  produces  liquefaction,  they  have  not  yielded 
diagnostic  liquefaction  parameters  that  are  easily  measured  in  situ. 

Currently,  in  situ  estimation  of  liquefaction  potential  is  accomplished 
using  procedures  based  on  the  standard  penetration  test,  cone  tip 
penetration  resistance,  and  normalized  shear  wave  velocity  (Finn  1995a,b). 
Such  procedures  involve  three  steps  (after  Finn): 

1)  Characterizing  the  dynamic  effects  of  the  earthquake. 

2)  Characterizing  the  in  situ  state  of  the  soil. 

3)  Application  of  a  criterion  for  the  incidence  of  soil 
liquefaction. 

Step  (1)  involves  equivalencing  of  anticipated  maximum  earthquake 
ground  motion  and  cyclic  laboratory  excitation  test  levels  and  durations. 
Step  (3)  involves  empirical  correlation  of  in  situ  test  probe  results  with  site 
specific  soil  materials  that  are  considered,  by  supplemental  criteria,  to  be  of 
high  liquefaction  potential.  Step  (2)  involves  development  of  definitive 
field  techniques  that  characterize  the  in  situ  state  of  the  soil  from  the 


reference  viewpoint  of  laboratory-generated,  engineering  criterion  of  soil 
liquefaction. 

This  report  considers  only  seismic  methods  of  imaging  subsurface  soil 
structure.  Primarily  focus  will  be  upon  Step(2),  and  secondarily  upon  Step 
(3).  The  specific  report  objectives  encompass  the  following  topics: 

1)  Evaluation  of  current  seismic  methods. 

2)  Postulation  of  potential  seismic  techniques. 

3)  Conceptual  definition  of  in  situ  field  tests. 

The  treatment  of  this  subject  begins  by  considering  soil  parameters  that 
are  relevant  to  the  liquefaction  problem.  Next,  since  the  process  of 
subsurface  imaging  requires  mathematical  models  of  soil  structure,  current 
models  for  lithified  and  unlithified,  saturated  and  unsaturated  media  are 
discussed.  A  comprehensive  treatment  of  seismic  tomography  follows, 
including  discussions  of  potential  capabilities  and  limitations.  Finally,  with 
reference  to  current  and  near-future  data  analysis  capability,  conceptual 
field  tests  are  presented . 
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2  Relevant  Soil  Parameters 


In  the  fluidization  process,  earthquake  wave  motion  imparts  pseudo- 
cyclic  loading  to  the  soil,  resulting  in  compaction.  Since  water  in  the  pores 
of  the  sediment  cannot  escape  quickly  enough  to  accommodate 
compaction,  an  increase  of  porewater  pressure  results,  leading  to  a 
reduction  of  intergrandular  stresses  between  soil  particles.  For  a  given 
intensity  of  earthquake  shaking,  liquefaction  potential  is  controlled  by  the 
tendency  for  compaction  and  capacity  for  drainage  (Finn  1995b). 

Cyclic  shearing  strain  is  recognized  as  the  prime  agent  for  generating 
volume  changes  (i.e.,  compaction)  in  saturated  cohesionless  soils.  Shear 
strain  depends  on  the  in  situ  shear  modulus,  which,  for  saturated 
cohesionless  soils,  is  strongly  dependent  upon  porosity.  Drainage,  on  the 
otherhand,  is  a  function  of  in  situ  permeability,  which  is  also  dependent 
upon  porosity.  Therefore,  porosity,  shear  modulus  and  permeability  are 
important  soil  parameters  for  the  evaluation  of  in  situ  liquefaction  potential 
(Finn  1995b). 

The  shear  modulus,  however,  cannot  generally  be  measured  directly. 
This  quantity  is  given  by  the  product  of  density  and  the  square  of  the  shear 
velocity.  Shear  velocity  can  be  readily  measured  in  situ.  If  in  situ  density 
was  also  determinable,  then  the  in  situ  shear  modulus  could  be  calculated. 
Density  is  also  dependent  upon  porosity  and  is  a  parameter  used  in  the 
calculation  of  ‘relative  density’.  Relative  density  is  a  liquefaction  indicator 
applied  to  the  evaluation  of  the  shearing  resistance  of  sands.  Thus,  density 
and  relative  density  are  also  liquefaction-relevant  soil  parameters. 

A  high  degree  of  fluid  saturation  is  required  for  liquefaction.  For  a 
partially  saturated  soil,  compaction  can  occur  without  significant  pore 
water  pressure  development.  Consequently,  the  in  situ  degree  of  saturation 
is  an  important  quantity  to  be  evaluated. 

Finally,  soil  fabric  and  structure  of  a  deposit  provide  clues  to  density, 
shearing  stiffness,  shearing  resistance,  and  porosity. 

In  summary,  the  following  soil  parameters  are  relevant  to  estimation  of 
liquefaction  potential: 

1)  density 

2)  relative  density 

3)  shear  modulus 


4)  porosity 

5)  permeability 

6)  saturation 

7)  soil  fabric  and  structure 


3  Seismic  Earth  Models 


Overview 


Geophysics  in  general,  and  seismology  in  particular,  is  an  interpretive 
science.  In  seismology,  earth  parameters  are  not  measured  directly,  but  are 
indirectly  infen'ed  through  comparison  of  recorded  waveform  attributes 
with  those  derived  synthetically  through  numerical  modeling. 
Consequently,  in  the  application  of  any  seismic  method  to  subsurface 
exploration,  one  must  of  necessity,  at  the  onset,  have  a  specific  earth  model 
in  mind,  compatible  with  a  chosen  process  for  numerical  synthesis  of 
seismic  data.  Due  to  the  complexity  of  numerical  modeling,  this  earth 
model  will  necessarily  be  simplistic;  inherently  capable,  at  best,  of 
reproducing  gross,  dominant  attributes  of  measured  seismic  waveforms. 
Realization  of  correspondence  between  field  measurements  and  numerical 
model  data  will  ultimately  depend  upon: 

1)  The  degree  of  inclusion  of  stress-strain  subtleties. 

2)  Significance  of  these  subtleties  relative  to  modification  of  the 
seismic  waveform. 

3)  V alidity  of  assumptions  of  isotorpy,  homogeneity  and  model 
dimensionality. 

4)  Accuracy  of  numerical  simulation  of  spatial  and  temporal 
source  characteristics. 


Models  For  Forward  Wavefield  Calculations 


Poroelastic  Models 

Biot  (1962a)  presents  a  unified  treatment  of  the  mechanics  of 
deformation  and  acoustic  propagation  in  porous  media.  He  treats  the  fluid- 
solid  rnedium  as  a  complex  physical-chemical  system  with  resultant 
relaxation  and  viscoelastic  properties  of  a  very  general  nature.  For  an 
isotropic,  elastic,  porous  medium,  the  principal  result  of  this  treatment  is 
the  field  equations  describing  coupled  energy  propagation  within  the 
porous  solid  structure  and  within  the  fluid.  These  equations  are  given  by 
Biot  as. 
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is  the  displacement  within  the  skeletal  solid, 
is  the  average  fluid  displacement, 
is  the  flow  of  the  fluid  relative  to  the  frame, 
is  the  dynamic  shear  modulus, 
is  the  density  of  the  fluid, 

is  theadiabatic  Lame  coefficient, 
is  the  porosity, 

is  the  complex  shear  modulus, 

is  the  specific  loss  of  the  skeletal  frame  for  shear, 

is  the  kinematic  viscosity  of  the  fluid, 

is  the  virtual  mass  of  the  skeletal  frame, 

is  the  total  mass  density, 

is  the  hydraulic  coefficient  of  permeability, 

are  Biot  elastic  coefficients  (Biot(1941)). 


Equation  (1)  is  the  result  of  the  application  of  the  principle  of  conservation 
of  momentum,  while  (2)  can  be  interpreted  as  expressing  the  dynamics  of 
relative  motion  of  the  fluid  in  a  frame  of  reference  moving  with  the  solid. 
These  equations  were  based  on  a  semiphenomenological  formulation  of  the 
equations  of  elasticity  for  a  porous  aggregate  by  Biot  (1941;  1956a,b,c; 
1962a,b),  Frenkel  (1944),  and  Gassmann  (1951).  The  hypothesized  form 
of  the  microscopic  constitutive  relationships,  which  resulted  in  the 
equations  of  motion,  was  confirmed  by  Burridge  and  Keller  (1981)  in  their 
studies  of  the  dynamic  equations  which  govern  the  behavior  of  the 
medium  on  the  microscopic  scale.  Moreover,  through  experimental 
observations,  Fiona  (1980)  confirmed  a  fundamental  prediction  of  Biot’s 
model  by  his  observation  of  a  second  compressional  wave  (referred  to  as 
the  ‘slow  wave’).  Stoll  and  Bryan  (1970)  and  Stoll  (1974)  extended  Biot’s 
theory  to  include  dissipation  due  to  the  frame  by  incorporating  complex 
moduli. 

Successful  applications  of  Biot’s  theory  to  various  fields,  notably  to 
seismic  exploration  (e.g.,  Geertsma  1957;  Geertsma  and  Smit  1961; 
Gardner,  Gardner,  and  Gregory  1974;  Domenico  1974;  Rosenbaum  1974;), 
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marine  acoustics  (Stoll  1974,  1977, 1980, 1985;  Yamamoto  1983;  Badiey 
and  Yamamoto  1985;  Turgut  and  Yamamoto  1988;  Yamamoto  and  Turgut 
1988;  Yamamoto  et  al.  1989),  and  attenuation  and  dispersion  of 
compressional  waves  (Dutta  and  Ode  1979,1983;  Stoll  and  Kan  1981; 
Ogushwitz  1985;  Geli ,  Bard,  and  Schmitt  1987;  Winkler  et  al.  1987; 
Rasolofosaon,  1987;  Schmitt,  Bouchon,  and  Bonnet  1988;  Yamamoto, 

Nye,  and  Kuru  1994,  1995)  have  demonstrated  the  utility  of  the  model  for  a 
variety  of  poroelastic  problems.  The  model  incorporates  eight  material 
parameters  (considering  77/ k  as  a  single  variable),  including  the 
liquefaction-sensitive  parameters:  density,  shear  modulus,  porosity,  and 
permeability.  Soil  fabric  and  stmcture  can  be  included  by  specification  of 
spatial  variation  of  material  coefficients.  Thus,  the  Biot  model  is  a  general 
formulation  of  small  amplitude  (linear  stress-strain  relation)  wave 
propagation  in  saturated  sediment  deposits  characterized  by  an  elastic, 
porous  frame.  In  the  forward  modeling  arena  (as  opposed  to 
interpretational  problems  of  waveform  inversion)  the  system  of  equations 
lends  itself  readily  to  waveform  calculations  by  means  of  finite  difference 
(Hassanzadeh  1991;  Zhu  and  McMechan  1991;  Dai,  Vafidis,  and 
Kanasewich  1995)  and  finite-element  methods  (Zienkiewicz  and  Shiomo 
1984). 


Biot-squirt  models 

In  Biot’s  model,  pore  fluid  is  forced  to  participate  in  a  solid’s  oscillatory 
motion  by  viscous  friction  and  inertial  coupling.  A  different  mechanism  of 
fluid  flow  during  seismic  and  acoustic  wave  propagation  is  associated  with 
the  squirting  of  pore  fluid  out  of  thin  soft  cracks  as  they  are  deformed  by 
passing  seismic  waves.  Mavko  and  Nur  (1979)  have  shown  that  the  squirt- 
flow  mechanism  results  in  much  higher  and  realistic  attenuation  values  in 
partially  saturated  rocks  than  those  predicted  by  Biot’s  mechanism.  Mavko 
and  Nur  (1975)  have  suggested  that  squirt  flow  may  occur  even  in  fully 
saturated  rocks  due  to  fluid  flowing  between  saturated  cracks  of  different 
orientation.  This  mechanism  has  been  shown  to  be  responsible  for  the 
measured  seismic  energy  losses  and  velocity  dispersion  in  sedimentary 
material  for  both  P-  and  S-waves  (e.g.,  Murphy,  Winkler,  and  Kleinberg 
1986;  Wang  and  Nur  1990).  Dvorkin  and  Nur  (1993)  and  Dvorkin, 
Mavko,  and  Nur  (1995)  proposed  a  combined  Biot-squirt  (BISQ)  model 
where  the  Biot  elastic  frame  is  replaced  with  a  viscoelastic  one.  They 
considered  a  macroscopically  homogeneous  rock  with  pore  space  that  has 
thin  compliant  (soft)  and  stiff  portions.  At  high  confining  pressure  the  thin 
compliant  pores  close  and  the  observed  velocity  dispersion  is  small  and  can 
be  adequately  described  by  Biot’s  theory  (Mavko  and  Jizba  1991). 
However,  at  small  confining  pressure  the  high  observed  velocity/frequency 
dispersion  and  attenuation  can  be  adequately  described  by  the  BISQ  model. 


Modeling  Summary 


A  BISQ  model  appears  to  be  appropriate  for  the  description  of  dynamic 
response  of  fully  or  partially  saturated  porous  media  to  low  amplitude 
seismic  wave  energy.  This  model  is  reducible,  by  certain  assumptions 
relative  to  the  independent  material  parameters  listed  above,  to  simpler  and 
standard  equations  of  motion  encountered  in  physics.  As  for  example, 
assuming  the  medium  to  be  elastic  and  non-porous,  (2)  above  vanishes, 
and  (1)  assumes  the  form, 


21^ 

j^j 


“ 

d 

1 

+  - 
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This  equation  can  be  recognized  as  the  isotropic,  heterogeneous  wave 
equation  for  a  perfectly  elastic  solid.  Other  equations  of  motion  can  be 
similarly  produced  by  appropriate  choices  of  material  parameter  behavior. 

I  conclude,  by  virtue  of  the  wide  application  of  Biot’s  theory  in  the 
literature  and  the  essential  absence  of  any  other  widely  accepted  model, 
that  the  BISQ  model  can  be  considered  the  starting  point  for  most  scientific 
inquiries  of  wave  motion  in  porous,  partially  or  fully  saturated,  earth 
materials.  The  equations  of  motion,  (1)  and  (2),  even  modified  by 
inclusion  of  squirt  flow,  are  amenable  to  numerical  solution  by  Finite 
difference  methods  (e.g.,  Dai,  Vafidis,  and  Kanasewich  1995). 
Consequently,  forward  modeling  of  stress  waves  in  porous  media  is  a 
matter  of  application  of  known  and  proven  numerical  methods,  and  is  not  a 
subject  requiring  extensive  research  and  development. 


4  Inversion  of  Seismic  Data 


Limited  Scope  Of  inversion  Methods 


In  the  discussion  of  seismic  inversion  techniques,  topics  of 
consideration  are  restrained  to  the  imaging  of  very  near-surface  materials 
typically  associated  with  earthen  dam  sites.  Complex  soil  structures,  with 
significant  variation  of  material  parameters  both  laterally  and  vertically,  are 
assumed.  Concepts  of  uniform  layered  media  are  not  applicable. 
Consequently,  imaging  and  migration  techniques  based  upon  layered  model 
concepts  are  not  considered. 

Seismic  inversion  methods  included  in  the  study  of  liquefaction 
potential  involved  direct  comparison  of  measured  waveform  attributes  with 
corresponding  waveforms  derived  synthetically  from  a  numerical  model. 
These  attributes  may  simply  be  traveltime  of  the  first  arrival,  or  first  arrival 
amplitude,  or  phase;  they  could  also  be  the  entire  pulse  waveforms,  or 
waveform  trains.  Regardless  of  character,  conformance  between  measured 
and  synthetic  data  attributes  generally  results  from  iterative  manipulation  of 
earth  model  parameters.  This  manipulation  is  exercised  in  a  regular  and 
controlled  way  by  the  minimization  of  some  objective  data  function;  the 
form  of  which  varies  in  accordance  with  the  particular  waveform  attribute 
in  question.  This  manipulation  process  is  called  tomography. 

Tomographic  algorithms  generally  fall  into  three  classes:  waveform, 
diffraction,  or  traveltime/attenuation  tomography.  For  waveform 
tomography,  synthetic  waveform  calculations,  based  on  some  current 
material  parameter  distribution,  must  be  made  by  finite  difference 
methods,  or  some  other  appropriate  modeling  procedure.  Waveform 
tomography  is  then  the  interactive  minimization  of  an  objective  function 
defined  in  terms  of  the  sum  of  time  integrals  of  differences  between 
observed  waveforms  and  synthetic  waveforms  calculated  by  a  forward 
modeling  procedure. 

For  diffraction  tomography,  synthetic  seismograms  are  generated  by 
finite  difference  methods,  or  some  other  appropriate  modeling  procedure, 
for  a  background  medium  void  of  anomalous  parameter  distributions.  The 
objective  function  is  then  the  difference  between  the  uniform  background 
velocity  and  the  anomalous  velocity  distribution  that  is  calculated  by 
inversion  of  the  scattered  field  data. 


In  traveltime  (or  attenuation)  tomography  the  forward  modeling 
problem  is  relatively  simple  and  direct,  involving  traveltime  (or 
attenuation)  calculations  through  a  grid  of  cells,  within  which  wave  speed 
(or  attenuation)  is  specified,  by  means  of  ray  or  wavefront  theory. 
Traveltime  (or  attenuation)  tomography  is  then  the  iterative  minimization 
of  an  objective  function,  defined  in  terms  of  the  sum  of  differences  between 
observed  traveltimes  (or  amplitudes)  and  traveltimes  (or  amplitudes) 
calculate  by  means  of  a  forward  model,  for  current  estimates  of  the  material 
parameters. 

For  a  comprehensive  presentation  of  tomographic  methods  the  reader  is 
referred  to  Nolet  (1987a),  Wu  and  Toksok  (1987),  and  Tien  and 
Inderwiesen  (1994). 


Waveform  Tomography 

Concept  Formulation 


I  begin  an  overview  of  tomography  by  considering  first  the  most 
general,  and  perhaps  the  most  potentially  powerful,  tomographic  process. 
Waveform  tomography  attempts  the  extraction  of  earth  parameters  from 
seismic  data  by  the  process  of  fitting  synthetic  wavetrain  elements  to 
observed  wavetrain  elements  within  a  measured  seismic  trace,  or  collection 
of  traces  (a  seismogram),  or  a  collection  of  seismograms.  Such  a  process, 
in  principal,  encompasses  the  possibility  of  determining  all  material 
parameters  of  the  model  chosen  as  the  best  representation  of  the  real  earth. 
Thus,  in  principal,  the  possibility  exists  to  simultaneously  determine 
compressional  and  shear  wave  velocities  and  density  as  a  function  of 
position  for  the  elastic  wave  model  (Equation  (3)),  or  all  eight  independent, 
spatially  variable,  material  constants  for  the  Biot  model  (Equations  (1)  and 
(2)).  I  stress  the  qualifier  in  principal’  in  that  the  quality  of  data,  relative 
significance  of  parameter  influence  on  the  shape  of  the  wavetrain,  and  the 
inaccuracy  of  inversion  procedures  limit  current  inversion  applications. 

Waveform  fitting  can  be  described  as  an  inverse  problem  based  upon 
methods  of  nonlinear  optimization.  In  general  the  optimization  problem  is 
one  of  minimizing  some  nonlinear  objective  function  of  the  model 
parameters,  which  is  generally  taken  to  be  the  difference  between  observed 
and  predicted  seismograms.  Inversion  is  then  reduced  to  the  question  of 
how  to  minimize  the  objective  function.  The  following  is  a  heuristic 
presentation  of  waveform  tomography. 

A  seismogram  is  acquired  as  output  of  discrete  sensors  distributed 
somewhere  in  space.  This  data  is  the  result  of  the  response  of  these  sensors 
to  incident  energy  generated  at  some  other  discrete  point  in  the  same  space. 
This  space  could  be  one-,  two-,  or  three-dimensional. 

Let  p  be  an  M-dimensional  vector  of  all  model  parameters  (velocity, 
density,  porosity,  permeability,  etc.);  the  i-th  synthetic  time  series 

calculated  by  means  of  model  p;  N,  the  total  number  of  seismic  time  series 
available;  T),  a  (sufficiently  large)  time  span;  si(t),  the  i-th  observed 


seismic  time  series;  and  q,  the  norm.  Following  Nolet  (1987b),  the 
objective  function  to  be  minimized  is  defined  by 

^(P)  =  -  S  J 1^  (P>  ,  (4) 

9  /=i  0 

where  R  is  a  desired  filtering  and/or  windowing  operator.  In  the  case  of 
traveltime  tomography  the  operator  R  would  be  unity  and  the  time  integral 
would  vanish.  Other  definitions  of  the  objective  function  may  be  used. 

The  reader  is  referred  to  Tarantola  (1987a). 

The  inversion  process  requires  certain  derivatives  with  respect  to  model 
parameters  pj.  First,  we  form  a  vector  whose  elements  are  given  by  the 
derivative  <3F(p)  /  Bp. , 

g(p)  s  VF(p)  =  X  - Ry,(0r''^^^ .  (5) 

1=1  0 

Next  we  form  the  Hessian  matrix,  H(p),  the  matrix  of  second  derivatives 
H(p)  =  Vg(p)  =  X  J{(^-  l)RVv^,[RV'/(P.O-R^,(Or'(RVv^,)'' 

i=l  0 

(6) 

+RVV  vr.[Rvr,.(p,0  -  Ry,.(0f 

Nonlinear  optimization  algorithms  work  in  an  iterative  way.  A  starting 
model  p  is  updated  with  a  correction  factor  Ap  to  give  a  new  model  p+ Ap. 
This  new  model  in  then  taken  as  a  starting  model  in  the  next  iteration. 
Suppose  that  at  some  stage  we  have  arrived  at  the  model  p.  With  a  simple 
Taylor  expansion  we  may  now  find  an  approximation  for  the  step  Ap  that 
should  bring  us  close  to  the  minimum  of  F(p),  where  g(p)  =  0: 

F(p  +  Ap)  =  F(p)  +  g(p)'^  Ap  +  ^  Ap'^H(p)  Ap .  (7) 

Differentiating  (7)  with  respect  to  Ap  gives: 

g(p  +  Ap)  -  g(p)  +  H(p) Ap .  (8) 

So  that  Ap  can  be  found  by  setting  g(p  +  Ap)  =  Oand  solving: 

H(p)Ap  =  -g(p).  (9) 

This  is  the  essence  of  waveform  tomography,  however,  it  is  noted  that 
the  above  presentation  is  given  heuristically  (with  much  liberality  relative 
to  mathematical  rigor).  This  was  done  for  the  purpose  of  clarity  of  process. 


For  a  more  rigorous  treatment  the  reader  is  once  again  referred  to  Nolet 
(1987b)  and  Tarantola  (1987). 


For  the  purpose  of  illustrating  how  the  vector  g(p)  and  the  Hessian 
matrix  H(p)are  calculated,  consider  the  heterogeneous  equations  of  motion 
for  a  solid  (non-porous)  media,  as  given  by  (3); 


i 


=  0. 


(3) 


In  this  equation  there  are  three  material  parameters:  //,  and  p,  which, 
in  general,  are  a  function  of  the  spatial  coordinates.  Consider  further  that 
we  approach  a  solution  to  this  equation,  and  the  partial  differential 
equations  to  follow,  by  means  of  finite  differences  (see  Kelly  (1976)). 

The  elements  of  g(p)  =  VF(p)  are  found  by  applying  the  operators 
d/dfi,  dldX^,  and  d!  dp  successively  to  Equation  (3),  resulting  in; 
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The  terms  in  the  R.H.S.  of  (10)  are  virtual  sources  generating, 
respectively,  'virtual  displacements': 

dy/j  _  du-  dyfj  _  du-  dy/j  _  du^ 

~^~~dp'  'd^~'^/  li^~'dp' 


(11) 


(since  de^ldp  =  d^y/j/dpdxj ,  etc.),  which  are  the  components  of  the  vector 

Vy/  in  (5)  and  (6).  By  first  solving  (3)  for  M,(x,t),  for  a  given  time  step 
within  the  finite  difference  iterative  time  loop,  the  virtual  source  terms  of 
(10)  can  be  readily  formed.  The  solutions  for  the  dependent  variables 


(virtual  displacements)  (11)  can  then  be  obtained  and  the  components  of 
g(p)  formed  for  the  same  time  step. 

To  calculate  the  elements  of  the  Hessian  matrix,  H(p)  ,  we  take 
derivatives  of  (10)  with  respect  to  all  three  parameters:  fi,  and  p.  For 
example,  one  derivative  produces  equations  of  motion  for  the  'virtual 
displacement'  d^xj/Jdpdu  =  d^ujdpdp: 
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The  terms  on  the  R.H.S.  of  (12)  are  virtual  sources  generating  virtual 
displacements:  / dpdu  =  d^u^l dpdp.  The  virtual  source  terms  of  (12) 

can  be  successively  formed,  (still  within  the  same  time  step)  from  the 
virtual  displacement  solutions  (1 1)  of  equation  (10).  Similar  equations 
would  exist,  in  this  example,  for  the  other  five  independent  components 
which  contribute  to  the  elements  of  the  symmetric  Hessian  matrix.  After 
all  ten  finite  difference  expressions  corresponding  to  Equations  (3), 
(lOa^c),  (12a ->f)  have  been  processed,  the  integrands  of  (4),  (5),  and  (6) 
are  formed  for  the  current  time  step  at  specific  positions  corresponding  to 
locations  of  receivers.  These  data  are  stored  and  the  loop  counter  is  then 
advanced  one  time  step.  This  looping  process  continues  until  a  set 
traveltime  Ti  is  reached  and  the  time  integrations  of  (4),  (5),  and  (6)  can  be 
completed. 

From  this  illustration  one  can  acquire  a  notion  of  how  waveform 
tomography  interrogates  the  independent  material  constants  in  reference  to 
their  individual  contribution  to  the  overall  waveform  under  scrutiny.  In 
this  finite  difference  example,  each  cell  of  the  numerical  model  would  be 
allowed  a  distinct  set  of  material  parameters  p,  ,  and  p  .  Because  of 
the  heterogeneous  formulation  of  the  equations  of  motion  (3),  continuity 
conditions  across  cell  boundaries  are  handled  implicitly,  allowing  for 
consideration  of  very  complex  models.  Only  external  boundaries  (free 
surface  and  artificial  model  bounds  at  the  bottom  and  sides)  are  handled 
explicitly.  Because  of  the  implicit  internal  boundary  conditions,  the  finite 
difference  method  is  a  common  method  of  choice  for  application  to  models 
with  complex  geometries. 

In  principle,  the  number  of  variable  material  parameters  (three  in  the 
previous  example)  is  not  limited.  A  tomographic  formulation  for  the  eight 
independent  material  parameters  of  the  Biot  model  can  be  similarly 
established.  However,  whether  a  stable,  meaningful  solution  of  the 
tomographic  process  for  the  full  Biot  model  is  achievable  or  not  is  very 
doubtful.  A  material  constant  hierarchy  exists  relative  to  significance  of 
their  effect  on  waveform  shape  and  character.  For  the  Biot  model  this 


ordering  has  not  yet  been  established.  Least  significant  parameters  would 
probably  be  masked  by  data  and  numerical  tomographic  process  errors. 


Application  State  of  Waveform  Tomography 

Full  waveform  tomography,  expressed  in  terms  of  nonlinear 
optimization,  is  not  limited  in  application  to  particular  data  acquisition 
geometries.  The  process,  in  principal,  can  be  applied  to  surface  reflection, 
VSP,  crosswell,  and  surface  wave  data.  In  the  literature,  however, 
applications  appear  only  for  the  inversion  of  global  surface  wave  data 
(Nolet  1987b)  and  for  inversion  of  multi-offset  seismic  reflection  data  (e.g., 
Bamberger  et  al.  1982;  Lines  and  Tritel  1984;  Tarantola  1984,  1986,  1987; 
Gauthier  Tarantola,  and  Virieux  1986;  McAulay  1985;  Kolb,  Collino,  and 
Lailly  1986;  Mora  1987,  1988;  Pan,  Phinney,  and  Odom  1988;  Pan  and 
Phinney  1989;  Helgesen  and  Kolb  1989;  Cao  et  al.  1990;  Pica,  Tarantola, 
and  Diet  1990;  Sen  and  Stoffa  1991;  Symes  and  Carazzone  1991;  Bunks  et 
al.  1995).  References  to  the  crosswell  problem  and  VSP  have  not  been 
found. 

For  waveform  tomography,  considerable  preprocessing  of  data  is 
required.  The  preprocessing  steps  are:  (1)  source  specification,  (2)  source 
radiation  correction,  (3)  attenuation  correction,  (4)  muting  and  windowing, 
(5)  wavelet  deconvolution,  and  (6)  two-and-a-half  dimensional  (2.5-D) 
corrections.  Preprocessing  procedures  are  given  by  Tura  et  al.  (1992). 

Most  iterative  full  waveform  inversion  methods  have  been  unsuccessful 
at  inverting  seismic  data  obtained  from  complicated  models.  The  primary 
difficulty  is  the  presence  of  numerous  local  minima  in  the  objective 
function  F(p).  The  presence  of  local  minima  at  all  scales  in  the  seismic 
inversion  problem  prevent  iterative  methods  of  inversion  from  attaining  a 
reasonable  degree  of  convergence  to  the  neighborhood  of  the  global 
minimum  (absolute  minimum  of  F(p)).  In  addition,  when  an  initial 
parameter  model  gives  rise  to  large  kinematic  errors  between  observed  and 
modeled  data,  a  perturbation  of  the  parameter  model  has  no  effect  on  the 
value  of  the  objective  function,  and  consequently  the  gradient,  g(p) ,  is  zero 
without  being  at  the  global  minimum.  The  solution  to  both  of  these 
problems  is  to  choose  an  initial  parameter  distribution  which  is  close  to  the 
global  minimum.  Obviously,  this  answer  is  not  the  answer,  since  the 
global  minimum  is  the  solution  that  we  seek.  But  the  problem  of 
convergence  to  the  global  minima  is  the  subject  of  intensive  current 
research  (e.g..  Bunk  1995). 

At  this  time,  waveform  tomography  is  more  of  an  art  form  than  a 
robust,  user  friendly  process.  Most  researchers  are  still  trying  to 
overcoming  the  local  minimum  problem  for  simply  velocity  models.  Mora 
(1987),  and  Sen  and  Stoffa  (1991)  have  extended  their  studies  to 
simultaneous  inversion  for  both  velocity  and  density,  while  Cao  et  al. 
(1990)  have  simultaneously  extracted  velocity  and  impedance.  But  it  is 
obvious  that  we  are  still  at  the  threshold  of  routine  waveform  tomographic 
applications,  and  certainly  not  close  to  consideration  of  multiparameter 
Biot  models.  While  the  process  machinery  is  essentially  in  place,  the  road 
map  for  routine  convergence  to  the  global  minimum  is  not. 


Diffraction  Tomography 

Concept  Formulation 

For  crosswell,  full  wave  imaging,  diffraction  tomography  appears  to  be 
the  method  of  choice.  This  method  uses  transmitted,  reflected,  and 
diffracted  events  to  provide  information  relative  to  material  parameter 
distribution  within  a  panel  bounded  by  two  boreholes.  The  method  also  has 
been  applied  to  VSP  data. 

The  basic  principle  of  diffraction  tomography  is  simple.  Consider  a 
wave  incident  on  an  object  in  a  homogeneous,  infinite  medium.  For 
simplicity,  consider  the  case  of  the  acoustic  wave  equation  with  constant 

density.  The  object  is  described  by  the  velocity  distribution  C(r),  where  r 
is  the  position  vector.  The  host  medium  has  a  velocity  Cq  .  The  wave 
equation  in  the  source-free  region  is 

VMr)  +  -^«(r)  =  0,  (13) 

C\r) 


where  «(r)  is  a  scalar  quantity  of  the  field  such  as  pressure,  co  is  the 
angular  frequency,  and  is  the  Laplacian  operator. 

Define  the  object  function  (9(r)  as 


0(r)  =  l- 


and  substitute  into  (13)  to  obtain 


(14) 


V^M(r)  +  k^u{r)  =  k^O{r)u{r) ,  (15) 

where  k  =  (o/Cq  is  the  wavenumber  of  the  field  in  the  host  medium.  Let 

M(r)  =  M°(r)  +  (7(r),  (16) 

where  M'’(r)  is  the  incident  wave  and  £/(r)  is  the  scattered  wave. 
Substituting  into  (15),  we  have 


V^U{r)  +  k^Uir)  =  eo{r)u{r). 


(17) 


By  using  the  free-space  Green’s  function  G(|r-  rj),  we  obtain 


(18) 


U{r)  =  -J  k^O{r'  )u{t'  )G(|r  -  r'|)rfr' , 

V 

where  the  integration  is  taken  over  the  volume  of  the  object.  Assuming  the 
object  is  a  weak  inhomogeneity,  the  Bom  approximation  (m  =  applies 
and  (18)  becomes 

t/(r)  =  -j  eO(r')G{\T  -  rj)«“(r')Jr’ .  (19) 

V 

Within  this  formulation  the  anomalous  velocity  variations,  as  expressed 
by  0(r),  form  a  basis  for  virtual  sources  impressed  on  an  otherwise 
homogeneous  medium  of  constant  velocity  Q.  These  virtual  sources  result 
in  a  scattered  field  U(t)  .  Diffraction  tomography  involves  the  inversion 

of  (19)  for  the  velocity  distribution  (9(r)  in  terms  of  U{t)  .  This  inversion 
is  accomplished  by  initially  assuming  that  the  incident  energy  is  a  plane 
wave,  and  that  the  receiver  is  sufficiently  removed  from  the  scattering 
objects  that  the  scattered  wave  may  be  considered  plane  also.  These 
assumptions,  together  with  the  Fraunhofer  approximation  for  the  Green’s 
function,  allows  the  scattered  plane  wave  in  a  given  direction  to  be  found 
from  the  3-D  Fourier  transform  of  the  object  function  (e.g.,  Devaney  1984). 
Theory  can  be  expanded  to  include  multiple  point  sources  (Harris  1987)  for 
reflection,  crosswell,  and  VS P  applications  (Wu  and  Toksbz  1987).  The 
problem  of  uniform  medium  assumption  has  been  addressed  by  Devaney 
and  Zhang  (1991),  Dickens  (1994),  Gelius  (1991,  1995a,b),  and  others. 


Application  State  of  Diffraction  Tomography 

Several  workers  have  used  synthetic  data  (Devaney  1982,  1984;  Wu 
and  Toksbz  1987;  Lo,  Duckwordi,  and  Toksbz  1990)  to  demonstrate  the 
ability  of  diffraction  tomography  to  produce  velocity  images,  with  a  spatial 
resolution  of  less  than  one  wavelength,  of  isolated,  weakly  scattering 
targets  embedded  in  a  constant- velocity  background.  Others  (Lo  et  al. 

1988;  Pratt  and  Worthington  1988)  have  used  diffraction  tomography  to 
image  low-contrast  scale  models,  and  the  algorithm  has  also  been  applied 
to  field  data  (Tura  et  al.  1992).  These  studies  have  used  the  filtered  back- 
projection  diffraction  tomography  algorithm  (Devaney  1982,  1984;  Wu  and 
Toksbz  1987)  which  is  applicable  only  for  the  case  of  weak  scatterers  in  a 
constant-velocity  medium.  In  general,  these  restrictions  make  it  impossible 
to  apply  the  traditional  filtered  back-propagation  diffracted  tomography 
algorithms  to  arbitrary  models  of  realistic  geological  complexity. 

Dickens  (1994)  applied  diffraction  tomography  to  the  problem  where 
the  velocity  structure  can  be  approximated  by  a  set  of  horizontal  layers  He 
shows  that  given  an  initial  layered  model  that  adequately  represents  the 
average  structure  of  the  subsurface,  and  an  accurate  calculation  of  the 
background  wavefield  propagating  through  this  model,  layered  diffraction 
tomography  can  be  used  to  successfully  image  complex,  geologically 
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realistic  models  to  which  the  traditional  filtered  back-propagation 
diffraction  topography  algorithm  is  inapplicable. 

A  generalized  diffraction  tomography  algorithm,  applicable  for  an 
arbitrary  background  structure,  has  been  studied  by  Gelius  et  al.  (1991), 
and  Gelius  (1995a,b).  The  method  is  valid  for  a  two-dimensional 
nonuniform  background  model  and  point  source  illumination  (i.e.,  a  2.5- 
dimensional  geometry).  This  method  is  regarded  as  a  generalization  of  the 
iterative  algorithm  of  Ladas  and  Devaney  (1991),  which  is  valid  only  for 
line  sources  and  two-dimensional  homogeneous  background  models. 

As  in  the  case  of  waveform  tomography,  extensive  preprocessing  of 
field  data  is  required.  And,  as  in  the  former  case,  only  first  arrival 
waveform  energies  are  generally  considered.  The  necessity  of  assumption 
of  weak  scattereres  in  a  constant- velocity  medium  eliminates  traditional 
filtered  back-projection  diffraction  algorithms  for  application  consideration 
to  the  liquefaction  problem.  On  the  other  hand,  the  diffraction  algorithms 
of  Gelius  (1995b)  look  promising. 

Categorically,  diffraction  tomography  is  based  on  the  acoustic  wave 
equation  and  has  focused  only  on  inversion  for  anomalous  velocity 
distributions.  There  are  inherent  difficulties  within  the  formulation  that 
would  be  presented  if  density  variations  were  also  allowed.  Thus,  density 
is  always  assumed  constant,  and  only  first  arrival  P-wave  energy 
waveforms  are  used  in  the  formulation  to  obtain  estimates  of  compressional 
wave  velocity  distributions. 


Traveltime/Attenuation  Tomography 

Concept  Formulation 

Traveltime/attenuation  tomography  begins  with  establishing  a  cellular 
model,  which  dimensionally  corresponds  to  a  section  of  the  real  earth  under 
investigation.  Each  cell  of  the  model  is  a  square  (or  triangle)  in  two- 
dimensions,  or  a  cube  (or  tetrahedron)  in  three-dimensions.  Associated 
with  each  cell  is  an  assigned  slowness,  sjj.  Provisions  must  then  be 
established  for  calculation  of  minimum  traveltimes  from  source  positions  to 
receiver  positions  for  waves  traversing  the  cellular  structure.  For  each 
source-receiver  pair  the  minimum  travel  path  must  be  superposed  upon  the 
cellular  structure,  and  the  path  length  segment  must  be  calculated  for  each 
cell  intercepted  by  the  ‘ray’.  These  path  length  segments  form  the  matrix 

A. 


A  common  approach  to  traveltime/attenuation  tomography  begins  with 
the  basic  equation 

N 

h  =  i  =  (20) 

j=i 
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where  r,.  is  the  traveltime  or  amplitude  associated  with  the  i-th  ray;  A.,  the 
pathlength  of  the  i-th  ray  in  the  j-th  cell;  and  Sj  is  the  slowness  or 
absorption  value  of  the  j-th  cell,  or  in  matrix  notation, 


t  =  As.  (21) 

When  (21)  is  solved  directly,  the  squared  solution  length  is  given  by  s’^s, 
and  constraints  involving  solution  length  seek  a  model  with  minimum 
overall  parameter  estimates.  This  is  not  desirable,  and  provides  motivation 
for  reformulating  the  problem  in  terms  of  a  background  model  and 
perturbation  values,  exactly  as  is  done  when  linearizing  a  nonlinear 
problem.  Then  the  basic  equation  to  solve  is 


At  =  AAs,  (22) 

where  As  =  s  -  ,  At  =  t  - 1^,  and  t,  are  vectors  of  theoretical  data  values 
associated  with  the  reference  model  .  Any  of  several  methods  might  be 

used  to  solve  (22)  for  As ,  including  the  direct  or  singular  value 
decomposition  (SVD)  computed  generalized  inverse  (Golub  and  Reinsch 
1970;  Lines  and  Treitel  1984;  Lines  and  Lafehr  1989),  back-projection 
methods  such  as  ART  (Tanabe  1971;  Dines  and  Lytle  1979)  or  SIRT 
(Gilbert  1972;  Ivansson  1986),  sparse  system  conjugate  gradient  methods 
(Hestenes  and  Stiefel  1952;  Paige  and  Saunders  1982;  Scales, 
Gersztenkom,  and  Treitel  1988),  and  composite  distribution  methods 
(Clippard,  Christensen,  and  Rechtien  1995).  For  most  tomography 

problems,  A  is  ill-conditioned,  and  direct  solution  results  in  unstable 
model  parameter  estimates.  To  stabilize  the  solution,  a  priori  information 
must  be  added.  This  information  in  usually  specified  in  the  form  of 
auxiliary  constraints  implemented  via  weighting  matrices.  The  objective 
function  to  be  minimized  is  then  given  by 

£'  =  (AAs-At)  W^(aAs- At)-l-As’'W„As,  (23) 


where  and  are  model  and  error  weighting  matrices,  respectively.  If 

we  define  and  s  DJD^,  then  an  penalty  function 

interpretation  of  E  is 


A 

As- 

D  Atl 

e 

LOn, 

i 

0  J 

(24) 


The  weighted  least-squares  parameter  estimates  may  be  written  as 


As  =  (a^DJD.A  +  d;;d„)’' A’'DjD,At. 


(25) 
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The  success  of  this  inversion  process,  yielding  corrections  As  to  the 

background  slowness  s„ ,  depends  upon  the  accuracy  of  the  elements  of  A. , 
which  in  turn  depend  upon  the  method  used  to  calculate  minimum 
traveltimes.  Traveltimes  may  be  obtained  by  direct  (numerical)  or  indirect 
(ray  tracing)  solution  of  the  eikonal  equation  (see  Cerveny  and  Hron  1980). 
The  only  field  variable  required  to  solve  the  eikonal  equation  is  the 
propagation  velocity,  but  direct  efficient  numerical  solution  of  the  eikonal 
equation  is  difficult  when  the  velocity  field  is  complicated.  Thus,  many 
traveltime  applications  involve  ray-tracing  methods.  Ray-tracing  equations 
are  ordinary  differential  equations  derived  by  applying  the  method  of 
characteristics  to  the  eikonal  equation  (Cerveny  and  Hron  1980).  The  ray 
equations  may  be  solved  with  shooting  methods,  which  pose  ray  tracing  as 
an  initial  value  problem,  or  with  bending  methods,  which  pose  ray  tracing 
as  a  two-point  boundary  value  problem. 

Alternately,  other  direct  numerical  solutions  of  the  eikonal  equation  are 
possible.  Vidale  (1990)  succeeded  with  a  finite  difference  scheme  that 
makes  a  plane-wave  approximation.  A  similar  approach  was  taken  by 
Podvin  and  Lecomte  (1991).  Moser  (1991)  presents  a  graph  theory 
approach;  Schneider  et  al.  (1992),  a  Fermat’s  principle  approach;  van  Trier 
and  Symes  (!991),  a  two-dimensional  upwind  finite  difference  approach; 
and  Schneider  (1995),  a  three-dimensional  upwind  finite  difference 
solution  of  the  eikonal  equation. 


Application  State  of  Traveltime/Attenuation  Tomography 

First  arrival  tomographic  algorithms,  far  too  many  to  list,  have  proven 
to  be  quite  useful  for  delineating  subsurface  structure.  Most  of  these 
algorithms  are  two-dimensional,  limiting  imaging  to  a  panel  defined 
(generally)  by  two  vertical  boreholes.  In  every  case  it  is  the  first  arrival 
within  a  seismic  trace  that  is  used  for  input  data.  Thus,  methods  to  date 
have  considered  only  compressional  wave  energy. 

High  resolution  inversion  requires  that  the  number  of  cells  (pixels) 
within  the  image  section  be  less  than  or  equal  to  the  number  of  independent 
traveltime  (or  amplitude)  measurements.  By  independence  it  is  implied, 
for  example,  that  two  closely  spaced  parallel  rays  traversing  the  same  cells 
would  only  count  as  a  single  measurement.  Thus,  the  available  number  of 
independent  traveltime  measurements  determines  the  size  of  the  image  cell, 
and,  consequently,  spatial  image  resolution.  However,  if  one  does  not 
require  high  resolution  of  image  parameters  (via  optimum  minimization  of 
the  objective  function),  but  rather  desires  a  more  detailed  relative 
distribution  image  (not  too  accurate  in  the  values  of  the  imaged  quantity), 
then  the  cell  size  can  be  made  smaller,  and  iterative  back-projection 
methods  such  as  ART  or  SIRT  can  be  used.  Examples  of  the  latter 
technique  is  given  by  Rechtien  and  Ballard  (1993). 


Hybrid  Tomography 

Biot  Theory  Applications 

Yamamoto,  Nye,  and  Kuru  (1994)  introduced  a  procedure  to  extract  the 
material  properties  of  unlithified  sediments  such  as  porosity,  permeability, 
and  shear  strength  directly  from  crosswell  compressional  wave  velocity 
images.  The  Biot  (1956a)  theory  and  an  empirical  relation  (Yamamoto  et 
al.  1989)  between  porosity  and  shear  modulus  were  used  to  extract  the 
porosity  images  and  the  shear  strength  images  from  the  velocity  images. 
This  procedure  was  successfully  applied  to  an  alluvial  sediment  strata  in 
Tokyo  Bay.  Yamamoto,  Nye,  and  Kuru  (1995),  extended  the  application 
of  this  procedure  to  the  case  of  mixed  lithified  and  unlithified  sediments 
comprising  a  Florida  limestone  aquifer.  They  adapted  the  unified  BISQ 
flow  theory  in  Dvorkin  and  Nur  (1993),  Dvorkin,  Nolen-Hoeksema,  and 
Nur  (1994),  and  Dvorkin,  Mavko,  and  Nur  (1995)  for  extraction  of 
permeability,  and  the  empirical  model  in  Han,  Nur,  and  Morgan  (1986)  for 
extraction  of  porosity  and  shear  strength. 


Application  state  of  hybrid  tomography 

The  results  of  Yamamoto,  Nye,  and  Kuru  (1994)  for  the  unlithified 
sediment  in  Tokyo  Bay  are  impressive  for  the  extraction  of  porosity  and 
shear  strength  images  from  compressional  wave  crosswell  data.  The 
extraction  of  a  permeability  image  was  less  successful  due  to  non- 
observance  of  clear  velocity-frequency  dispersion  within  the  experimental 
frequency  range  (1  kHz  -  10  kHz).  These  results  were  based  on  the 
adoption  of  the  Biot  low-frequency  theory  as  an  adequate  unlithified 
sediment  (alluvial  clays  and  sands)  model. 

The  results  of  Yamamoto,  Nye,  and  Kuru  (1995)  were  for  a  mixed 
lithified  and  unlithified  sediment.  In  this  experiment,  porosity  and  shear 
velocity  estimations  were  based  upon  a  specific  empirical  model  of  a 
lithified  sediment,  as  given  by  Han,  Nur,  and  Morgan  (1986).  This  model 
first  relates  porosity  to  compressional  wave  speed  (as  given  by  velocity 
imaging  of  first  arrival  field  data)  and  a  measure  of  clay  content  (measured 
or  assumed).  Secondly,  shear  wave  velocity  is  expressed  in  terms  of 
porosity  (as  calculated  in  the  first  step)  and  clay  content.  Shear  strength  is 
then  calculated  by  using  assumed,  or  measured,  values  of  grain  and  pore 
fluid  densities.  These  calculations  can  be  performed  for  low-frequency 
compressional  wavefield  data  providing  a  low-frequency  empirical  lithified 
sediment  model  is  established. 

Permeability  was  estimated  by  fitting  the  velocity  -frequency  dispersion 
data  of  the  Biot  (1956)  model  or  the  BISQ  model  (Dvorkin,  Nolen- 
Hoeksema,  and  Nur  1994)  with  corresponding  field  data  extractions. 
Experimentally,  Yamamoto,  Nye,  and  Nur  (1994,  1995)  employed  an 
acoustic  pulse  with  its  energy  concentrated  near  a  certain  frequency  by 
using  a  narrow-band  pseudo  random  binary  sequence  code  as  a  source 
signal.  Crosswell  measurements  through  the  same  cross-section  were 
repeated  many  times  using  different  frequencies  ranging  from  200  to  5000 


Hz.  The  inhomogeneous  structure  of  rock  permeability  was  indicated  by 
images  obtained  by  taking  differences  of  velocity  images  acquired  at 
different  carrier  frequencies. 


Non-Tomographic  Seismic  Methods 


Empirical  methods  have  been  developed  to  evaluate  liquefaction 
resistance  directly  from  shear  wave  velocity  (Bierschwale  and  Stokoe 
1984).  Conventional  crosswell  (non-imaging)  and  downhole  shear  wave 
methods  have  been  employed  (Stokoe  and  Hoar  1978;  Woods  1978;  Sirles 
1987)  to  obtain  depth  estimates  of  shear  wave  velocity.  Nazarian  and 
Stokoe  (1984)  developed  the  Spectral  Analysis  of  Surface  Waves  (SASW) 
procedure,  where  the  distribution  of  shear  wave  velocity  with  depth  is 
determined  by  phase  velocity  analysis  of  Rayleigh  waves.  Sirles  (1987) 
investigated  liquefaction  in  terms  of  amplitude  attenuation  of  shear  waves 
in  a  crosswell  experiment.  These  methods  are  all,  essentially,  one¬ 
dimensional  in  that  they  give  variation  of  shear  waves  with  depth  only 
along  a  vertical  line  arbitrarily  positioned  at  the  midpoint  of  the  two 
measurement  points.  Basically,  the  fundamental  assumption  common  to  all 
these  methods  is  that  the  local  earth  is  composed  of  uniform  horizontal 
layers.  Some  semblance  of  lateral  variation  of  shear  velocity  can  be 
obtained  by  shifting  the  experimental  arrangement  (whether  boreholes  are 
surface  instrumentation)  laterally  and  collecting  a  suite  of  data  along  a 
horizontal  profile.  But  these  methods  are  not,  in  any  true  sense,  two- 
dimensional  as  are  the  tomographic  methods  previously  discussed. 
Nevertheless,  these  methods  have  produced  shear  wave  velocity  estimates 
useful  to  the  liquefaction  problem. 

Compressional  and  shear  wave  refraction  methods  are  also  occasionally 
applied.  However,  the  hidden  layer  problem  (low  velocity  layer  beneath  a 
high  velocity  layer)  presents  a  very  real  and  commonly  encountered  pitfall 
that  often  renders  refraction  profile  data  meaningless.  Worse  yet,  the 
presence  of  this  condition  is  difficult  to  detect  and  usually  must  be 
determined  by  some  other  means. 

Since  these  conventional  methods  are  familiar  to  investigators  in  the 
field  of  liquefaction,  no  further  discussion  of  them  will  be  made.  The 
interested  reader  can  refer  to  the  references  cited. 


5  Conceptual  Field  Tests 


Potential  Seismic  Methods  for  Liquefaction  Studies 


Shear  waves  are  the  dominant  energy  fields  used  for  the  assessment  of 
liquefaction  potential.  Shear  wave  velocity,  by  virtue  of  its  relation  to  the 
shear  modulus,  is  an  indicator  that  is  directly  linked  to  the  development  of 
pore  water  pressure  in  soils  under  cyclic  loading.  However,  since  shear 
wave  energy  is  never  the  first  arrival  in  a  seismic  trace,  even  though  shear 
sources  and  cross-polarized  acquisition  procedures  are  employed,  it  is 
extremely  difficult  to  achieve  high  accuracy  in  measurement  of  arrival  time 
of  shear  waves.  The  error  in  arrival  time  is  considerably  greater  than  that 
for  first  arrival  compressional  energy,  and,  consequently,  substantial  error 
in  shear  wave  velocity  results. 

To  my  knowledge,  traveltime  tomography  has  never  been  applied  (in 
the  classical  sense)  to  the  problem  of  direct  extraction  of  shear  velocity 
distributions.  But  it  is  possible,  in  principal,  to  directly  extract  shear  wave 
information  by  use  of  waveform  tomography.  So  we  are  faced  with  the 
dilemma  that  procedural  problems  for  the  latter  must  still  be  worked  out  to 
make  the  method  useful,  and  the  former  method  works  well  for 
compressional  waves,  but  not  shear.  Thus,  at  this  time,  it  seems  that  we  are 
left  with  conventional  methods,  as  discussed  in  the  previous  section,  and 
compressional  wave  traveltime/attenuation  tomography  to  investigate  the 
liquefaction  problem. 

Traveltime/attenuation  tomography  for  first-arrival  compressional  wave 
energy  is  useful  for  the  study  of  liquefaction  for  several  reasons.  First, 
tomographic  compressional  wave  velocity  images  yield  information 
relative  to  overall  fabric  and  structure  of  the  medium.  Secondly, 
attenuation  tomography  provides  spatial  distribution  of  overall  (all 
frequencies)  material  attenuation.  Frequency  dependent  attenuation,  or, 
equivalently,  the  quality  factor  Q,  can  possibly  be  obtained  by  imaging 
narrow-band  filtered  first  arrival  amplitudes.  Thirdly,  shear  wave  velocity 
and  porosity  distributions  can  be  derived  from  compressional  wave 
traveltime/attenuation  tomographic  data  by  application  of  the  empirical 
method  of  Yamamoto,  Nye,  and  Kuru  (1994).  Depending  on  signal 
bandwidth,  permeability  estimates  might  also  be  obtained  from 
compressional  wave  data  by  means  of  the  empirical  method. 


Compressional  wave  tomography  can  be  used  in  conjunction  with 
conventional  procedures.  Spectral  analysis  of  surface  waves  (SASW)  can 
be  conducted  independent  of  any  tomographic  data  gathering  activity.  This 
method  results  in  estimation  of  shear  wave  velocity  with  depth.  In 
addition,  crosswell  and  downhole  shear  wave  studies  can  be  made  using  the 
same  boreholes  employed  for  tomographic  studies.  SASW,  conventional 
crosswell  and  downhole  shear  wave  studies  are,  essentially,  one¬ 
dimensional  in  that  velocity  is  assumed  laterally  constant,  but  variable 
with  depth.  Thus,  tomography  is  here  proposed  as  a  supplemental 
procedure  to  conventional  one-dimensional  methods. 

Although  currently  constrained  by  inadequacy  of  data  processing 
procedures,  waveform  tomography  remains  the  imaging  goal.  The  advent 
of  adequacy  is,  perhaps,  near  future,  and  it  seems  reasonable  to  execute 
appropriate  data  gathering  procedures  that  would  render  data  suitable  for 
future  analysis  by  advanced  methods.  Consequently,  in  addition  to 
conventional  SASW,  crosswell,  and  downhole  shear  wave  studies,  I 
propose  the  following  data  acquisition  procedures: 

1)  Combined  crosswell  &  borehole-to-surface  (CCBS),  full 
waveform  data  gathering,  using  compressional  wave  sources  and 
hydrophones. 

2)  Combined  crosswell  &  borehole-to-surface  (CCBS),  full 
waveform  data  gathering,  using  clamped  shear  wave  sources  and 
clamped  three-component  seismometers. 

3)  Full  waveform  Rayleigh  wave  (FRW)  study  using  a  vertical 
surface  impact  source,  eighteen  three-component  surface 
seismometers,  and  two  three-component  clamped  borehole 
seismometers. 

4)  Full  waveform  Love  wave  (FLW)  study  using  a  horizontal  surface 
impact  source,  eighteen  three-component  surface  seismometers, 
and  two  three-component  clamped  borehole  seismometers. 


Experimental  Geometry  -  CCBS 


Field  Configuration 

Figure  1  shows  the  proposed  experimental  arrangement.  Assume  a 
zone  of  anomalous  wave  velocity,  density,  porosity,  permeability,  etc. 
embedded  in  a  typical  soil  environment.  We  make  no  assumption  of 
material  distributions,  water  saturation,  water  table  depth,  etc.. 

A  system  of  four  boreholes  is  desirable,  although  two  would  suffice. 
Each  borehole  consists  of  a  4  inch-diameter  (minimum)  PVC  pipe.  This 
pipe  must  be  appropriately  grouted  to  the  host  medium.  The  depth  of  each 
borehole  will  be  approximately  three  times  the  anticipated  depth  of  the 
anticipated  liquefiable  zone  (3D).  The  horizontal  spacing  of  these  holes 


will  be  approximately  one-third  the  anticipated  depth  of  the  zone  (D/3).  In 
the  event  of  having  only  two  boreholes,  the  desired  borehole  spacing  would 
be  (D).  These  boreholes  must  be  sealed  in  order  to  enable  the 


Figure  1.  Common  source  acquisition  configuration  for  tomography. 

maintenance  of  a  full  column  of  water  in  each  PVC  pipe.  In  addition, 
shallow  holes,  approximately  two  feet  in  depth,  will  be  emplaced  along  the 
profile.  These  holes  will  be  lined  with  plastic  to  enable  the  maintenance  of 
a  shallow  water  column. 

Two  separate  data  sets  will  be  gathered.  The  first  set,  as  illustrated  in 
Figure  1 ,  will  employ  broad  band  hydrophones  as  the  sensing  elements  and 
a  high-frequency  sparker  system  as  the  energy  source.  This  system  is 
designated  as  the  Borehole  Imaging  and  Tomographic  System  (BITS),  a 
product  of  past  Waterways  Experiment  Station  (\^^S)  tomographic 
research.  A  description  of  BITS  can  be  found  in  Rechtien,  Hambacker,  and 
Ballard  (1993).  The  second  data  set  will  be  gathered  by  the  use  of  wall- 
clamped  three-component  geophones,  surface  three-component  geophones, 
and  a  wall-clamped  vertical  polarizing  shear  wave  source. 

For  the  collection  of  data  suitable  for  waveform  tomography  (having  in 
mind  future  data  processing  capability),  it  is  necessary  that  all  receiving 
instrumentation  have  similar  dynamic  response.  Consequently,  it  is  not 
permissible  to  mix  instrumentation;  i.e.,  borehole  hydrophones  and  surface 
geophones. 

During  acquisition  of  both  sets  of  data,  source  locations  will  be 
established  in  one  of  the  outer  boreholes  and  receiver  locations  established 
in  the  other  outer  borehole  and  along  the  ground  surface  between  these  two 
boreholes.  As  indicated  in  Figure  1,  seismic  energy  generated  at  a  given 
source  position  will  be  recorded  (non-simultaneously)  at  selected  receiver 


positions.  The  acquisition  procedure  employed  is  designed  to  acquire 
maximum  density  of  ray  coverage  in  the  upper  two-thirds  of  the  vertical 
panel  defined  by  the  outer  two  boreholes. 

The  inner  two  boreholes,  separated  by  a  spacing  of  D/3,  shall  be  used 
for  the  collection  of  conventional  crosswell  and  downhole  shear  wave  data 
for  comparison  with  tomographic  shear  velocity  images. 


BITS  Instrumentation 

The  BITS  system  records  data  within  the  frequency  bandwidth  of  10  Hz 
to  2  kHz  (an  anti-alias  filter  is  set  at  2  kHz).  The  signal  bandwidth  dictates 
the  spacing  of  instrumentation  on  the  surface  and  within  the  borehole. 
Within  the  borehole,  hydrophones  will  sense  tube  waves  superposed  upon 


Figure  2.  Constant  vertical  offset  data  gathering  configuration.  A  single 
hydrophone  and  a  single  source,  vertically  offset  by  an  angle  0,  are 
simultaneously  raised  a  distance  Az  between  data  acquisition  cycles. 

Fixed  surface  hydrophones  are  sampled  for  every  source  position. 

desired  crosswell  body  wave  transmissions.  Tube  waves  are  generated 
primarily  by  the  very  process  of  source  wavefield  incidence  upon  the 
borehole.  A  second  tube  wave  effect  occurs  in  the  data  as  a  result  of  tube 
wave  conversion  within  the  source  borehole,  with  the  converted  waves 
acting  as  secondary  sources.  By  employing  a  constant  vertical  offset  data 
collection  procedure,  where  both  source  and  receiver,  vertically  separated 
by  a  given,  fixed  distance,  are  simultaneously  raised  to  shallower  positions 
between  recordings  (see  Figure  2),  both  tube  wave  effects  are  manifested 
within  the  constant  vertical  offset  trace  record  as  events  falling  along  a  line 
with  inverse  slope  approximately  equal  to  the  tube  wave  velocity.  These 
effect  can  be  readily  removed  from  the  data  by  used  of  frequency- 


wavenumber  velocity  filters  provided  the  data  is  not  spatially  aliased. 
Spatial  aliasing  occurs  when, 

A  K 

Az<^,  (26) 


where  Az  is  the  sensor  spacing,  V,  is  the  tube  wave  velocity,  and  is 
the  maximum  frequency  with  non-trivial  signal  strength.  The  tube  wave 
velocity  is  in  the  neighborhood  of  4000  fps,  given  a  water-filled  borehole. 
Thus  for  borehole  hydrophones,  and  a  maximum  frequency  of  2  kHz,  both 
source  and  hydrophone  station  spacing  should  be  0.5  ft  or  less. 

A  third  tube  wave  effect  to  consider  is  that  the  tube  wave  velocity 
might  exceed  the  compressional  wave  velocity  in  the  cohesionless  soil. 
While  tube  wave  energy  falls  off  rapidly  with  distance  from  the  center  of 
the  tube,  it  would  nevertheless,  in  reality,  be  a  secondary  source  moving  at 
‘supersonic’  speeds  through  the  medium.  The  resulting  shock  wave 
process  is  somewhat  similar  to  igniting  a  vertical  string  of  primacord 
suspended  within  a  borehole,  or  to  an  airplane  plowing  through  the 
atmosphere.  Consequently,  it  is  necessary  to  use  mechanical  tube  wave 
attenuators  both  above  and  below  the  source  and  above  and  below  the 
receiver.  These  attenuators  are  commercially  available. 

A  second  reason  to  employ  tube  wave  attenuators  is  that  tube  wave 
energy  encountering  the  free  surface  is  not  only  partially  reflected,  but 
significant  tube  wave  energy  is  converted  to  body  and  surface  waves. 

Since  tube  wave  energy  will  arrive  at  the  surface  in  advance  of  transmitted 
body  waves,  converted  energy  might  be  recorded  at  surface  hydrophone 
positions  as  first  arrivals.  Consequently,  the  desired  transmitted  events  will 
be  masked  by  these  converted  waves. 

A  final  consideration  concerns  estimation  of  source  signature.  For  the 
execution  of  waveform  tomographic  data  processing  the  characteristics  of 
the  source  must  be  precisely  known.  Waveform  inversion  requires  model 
waveform  data  to  be  numerically  generated  and  compared  with  field  data 
(specified  within  the  objective  function).  Thus,  modeled  data  must  contain 
not  only  the  complex  spectral  content  of  the  source,  but  also  the  spatial 
radiation  pattern  effects  of  source  and  receiver  borehole  coupling.  From 
the  viewpoint  of  modeling,  the  boreholes  don’t  exist;  only  point  sources 
and  point  receivers  buried  within  the  medium. 

Source  characteristics  can  generally  be  gleaned  from  the  set  of  first 
arrival  waveforms.  These  waveforms  include  the  effects  of  source  spectral 
content,  source  radiation  and  receiver  reception  patterns,  and  signal 
modification  along  the  transmission  path.  The  source  spectral  content  can 
be  continually  monitored  by  a  spatially  fixed  hydrophone  mounted  within, 
and  near  the  top  of,  the  source  borehole.  Consequently,  source  spectral 
variations  can  be  corrected  from  these  monitored  signaJs.  Sets  of  constant 
vertical  offset  crosswell  data  yield  information  relative  to  source  radiation 
and  receiver  reception  at  various  ray  angles.  Long  wavelength  filtering  of 
these  data  will  remove  effects  of  spatially  small  material  inhomogenities  in 
the  data,  thus  bringing  clarity  to  source  radiation  and  receiver  reception 
geometries.  Finally,  attenuation  tomography  clearly  shows  borehole 
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coupling  problems  due  to  inadequate  grouting  between  the  PVC  pipe  and 
host  medium.  Amplitude  and  traveltime  corrections  to  affected  traces  in 
the  total  data  set  will  remove  these  undesirable  tomographic  anomalies. 


Shear  Instrumentation 

For  the  emplacement  of  wall-clamped  instrumentation,  the  borehole  can 
be  air-filled.  While  tube  waves  still  exist  in  and  around  an  air-filled 
borehole,  the  magnitude  of  such  energy  is  not  significant.  Acoustic 
mufflers  above  and  below  the  instrumentation  will  suffice  to  remove  any 
tube  wave  problem.  However,  spatial  filtering  of  the  data  may  still  be 
required  in  order  to  arrive  at  definition  of  source  radiation  and  receiver 
reception  patterns  (as  discussed  in  the  previous  paragraph  on  BITS 
instrumentation).  Therefore,  precautions  must  be  taken  to  insure  that  the 
data  is  not  spatially  aliased. 

In  consideration  of  spatial  aliasing  in  the  previous  section,  the  velocity 
specified  in  Equation  (26)  was  referred  to  as  the  tube  wave  velocity.  When 
viewing  a  multi-trace  seismogram,  a  linear  alignment  of  energy  can  often 
be  observed.  For  water-filled  borehole  data,  the  obvious  energy  alignments 
are  tube  waves.  The  inverse  slopes  of  these  alignments  yield  the  tube  wave 
velocity.  For  air-filled  boreholes  tube  waves  should  not  be  observable.  If 
they  exist  at  all  their  velocity  would  be  close  to  1100  fps.  Thus  for  a 
maximum  signal  frequency  of,  say,  200  Hz,  Equation  (26)  would  require 
instrumentation  spacing  of  about  1.3  ft.  Although  it  is  doubtful  that  such 
alignments  will  be  observable  in  air-filled  borehole  data,  it  must  be  recalled 
that  instrumentation  spacing  also  determines  the  pixel  size  of  the 
tomographic  image.  Consequently,  a  spacing  of  1  to  1.5  ft  would  be  most 
desirable  for  shear  wave  data  gathering.  Previous  comments  concerning 
source  signal  estimation  for  hydrophones  is  equally  applicable  for  clamped 
instrumentation. 


Experimental  Geometry  -  FRW  &  FLW 


Field  Configuration 

Employing  a  48  channel  recording  system,  18  stations  will  be 
established  along  the  dam,  with  the  borehole  configuration  near  the  center 
of  the  spread.  Station  spacing  will  be  in  the  neighborhood  of  5  feet.  A 
three-component  seismometer  will  be  installed  at  each  station. 

Simultaneously  employing  an  additional  6  channel  recording  system, 
two  three-component  clamped  geophones  will  be  installed  at  depths  of  D 
and  2D  within  a  single  borehole. 

For  the  Rayleigh  wave  study  (FRW),  energy  will  be  supplied 
alternately  at  both  ends  of  the  line  by  means  of  a  vertical  hammer  blow 
upon  a  plate.  For  the  Love  wave  study  (FLW),  energy  will  be  supplied  by 
alternately  striking  both  ends  of  a  vertically  loaded  wooden  beam  with  steel 


end  caps.  The  beam  will  be  oriented  with  long  axis  horizontal  and 
perpendicular  to  the  instrumentation  line. 

In  the  process  of  waveform  tomography  the  modeled  wave  field  is 
generated.  Such  wave  fields  include  waveforms  at  depth.  Consequently, 
having  borehole  measurements  of  the  shallower  roots  of  both  Rayleigh  and 
Love  waves  provides  a  means  of  assessing  the  accuracy  of  the  tomographic 
process. 


Conventional  Seismic  Tests 


The  inner  two  boreholes  of  Figure  1  are  for  conventional  downhole  and 
crosswell  testing  (Butler  and  Curro,  1981).  Crosswell  testing  involves 
horizontal  transmission  paths  only,  and  employs  clamped,  vertical 
polarized  shear  sources  and  three-component  wall-clamped  geophones. 
Downhole  testing  involves  surface  SH  cross-polarized  sources  and  three- 
component  wall-clamped  geophones. 

The  primary  reason  for  the  additional  two  boreholes  is  that  tomographic 
images  are  least  accurate  near  instrumentation  boreholes.  The  reason  for 
lack  of  accuracy  is  the  convergence  of  all  source-receiver  raypaths  to  single 
points  at  the  borehole  rather  then  having  rays  more  uniformly  spread  out 
within  the  pixel,  as  they  are  in  regions  more  central  to  the  image.  Thus 
image  data  is,  supposedly,  more  accurate  in  the  central  regions  of  the 
tomograph,  and  probably  more  correlatable  to  convention^  shear  wave 
estimates.  However,  this  is  a  minor  point,  and  confinement  to  two 
boreholes,  instead  of  four  is  satisfactory  provided  the  borehole  spacing  D  is 
small  enough  to  avoid  lack  of  signal  strength. 

In  the  event  that  the  two  inner  boreholes  are  not  available,  no  additional 
data  gathering  is  required  for  conventional  crosswell  and  downhole 
procedures.  Conventional  crosswell  and  downhole  data  is  included  within 
the  tomographic  data  sets. 

No  additional  data  gathering  is  required  for  the  SASW  procedure. 
Spectral  analysis  of  surface  waves  can  be  conducted  between  any  two 
stations  of  the  FRW  data  set.  Thus,  numerous  pairs  of  surface 
measurements,  with  a  variety  of  station  spacing,  will  be  available  for 
spectral  analysis. 


6  Summary 


The  BISQ  theory  is  acceptable  as  a  theoretical  dynamic  model  of 
unlithified,  or  partially  lithified,  saturated,  or  partially  saturated,  soil  when 
subjected  to  low  amplitude  seismic  wave  energy.  Finite  difference 
numerical  modeling  methods  exists  for  calculation  of  synthetic  wave  fields 
within  a  heterogeneous  Biot  solid.  Extension  to  the  BISQ  solid  is 
straightforward.  Compressional,  shear,  and  fluid  wave  components  are 
incorporated  within  the  formulation.  Existing  computer  algorithms  can  be 
expanded  for  three-dimensional  geometries. 

Full-waveform  tomographic  methods  are  potentially  applicable  to  Biot 
and  BISQ  solids.  Potential  exists  for  simultaneous  inversion  of  multiple 
model  parameters.  But  problems  concerning  convergence  of  the  objective 
function  to  the  global  minimum  have  not  yet  been  fully  resolved. 
Consequently,  full-waveform  tomography,  while  certainly  the  most 
desirable  tomographic  application,  is  not  yet  on-line  as  a  routine  inversion 
method. 

Diffraction  tomography  is  based  on  the  acoustic  wave  equation  and  is 
focused  only  on  inversion  for  anomalous  velocity  distributions.  Inherent 
difficulties  within  the  formulation  would  be  presented  if  density  variations 
were  also  allowed.  The  method  is  applicable  only  to  weak  scatterers  in  a 
medium  of  constant  background  velocity.  Some  movement  of  the  theory 
has  been  made  toward  application  to  an  arbitrary  background  structure.  In 
its  present  form  if  offers  no  advantage  over  traveltime  tomography,  and, 
henceforth,  should  not  be  seriously  considered  as  a  potential  method 
applicable  to  the  liquefaction  problem. 

Travel  time/attenuation  tomography  is  fully  developed.  Inversion  of 
first-arrival  traveltime  data  yields  spatial  estimates  of  a  quantity  called 
‘velocity’.  This  quantity  can  be  correlated  with  the  spatially  variable 
material  velocity  only  if  realistic  transmission  paths  are  employed  within 
the  forward  model  calculation  of  the  inversion  process.  The  use  of  straight 
raypath  tomographic  algorithms,  for  example,  can  never  yield  correct 
velocity  estimates  within  a  spatially  limited  low  velocity  zone.  Curved 
raypath  and  wavefront  algorithms  yield  better,  but  not  totally  accurate, 
results.  Inversion  of  first-arrival  amplitudes  yields  estimates  of  a  spatially 
distributed  quantity  call  ‘attenuation’.  This  quantity,  however,  has  no 
direct  correlation  to  fundamental  model  parameters  (i.e.,  frame  or  fluid 


viscosity,  etc.).  Attenuation  images  only  reveal  zones  of  overall,  broad¬ 
band  energy  diminution. 

The  hybrid  tomographic  method  of  Yamamoto,  Nye,  and  Kuru  (1994) 
define  porosity  and  shear  velocity  as  image  parameters.  Estimates  of  these 
parameters  are  derived  from  compressional  wave  crosswell  data.  These 
estimates  are  based  upon  empirical  relationships  derived  from  Biot’s  low- 
frequency  theory.  While  their  results  are  interesting,  correlation  with 
independent  estimates  of  porosity  and  shear  velocity  must  be  made  for 
validation. 

It  is  taken  for  granted  that  the  exhaustive  seismic  response  of  a 
liquefiable,  cohesionless  soil  can  be  fully  captured  within  appropriate  field 
measurements  of  seismic  waves.  Instrumentation  and  field  procedures  are 
sufficiently  developed  for  accurate,  comprehensive  acquisition  of  seismic 
data.  The  potential  for  data  quality  far  exceeds  current  ability  to  extract 
liquefiable-relevant  material  parameters  from  recorded  seismic  waveforms. 
Nevertheless,  it  is  expedient  to  collect  data  in  anticipation  of  resolution  of 
current  inversion  limitations. 

Currently,  direct  extraction  of  material  parameters  from  a  tomographic 
seismic  data  set  is  limited  to  compressional  wave  velocity  and  attenuation. 
Shear  wave  velocity  may  also  possibly  be  directly  extracted,  given 
sufficient  angularity  of  borehole  shear  source  radiation  strength  and  ability 
to  separate  out  compressional  energy  arrivals.  Empirical  procedures  allow 
the  expression  of  shear  wave  strength  and  porosity  from  compressional 
wave  data.  Permeability  estimates  may  possibly  be  empirically  derived, 
given  sufficient  data  bandwidth.  Compressional  wave  tomographs  also 
yield  information  relative  to  soil  fabric  and  structure. 

Surface  wave  studies  will  generate  data  suitable  for  future  waveform 
tomography.  Currently,  however,  conventional  SASW  techniques  can  be 
applied  to  both  Rayleigh  and  Love  wave  data  sets  (a  SASW  procedure  for 
Love  waves  is  probably  available  somewhere). 

Conventional  crosswell  and  downhole  shear  wave  tests  are  well 
established. 


7  Conclusion 


An  existing,  limited  data  set  suggest  that  shear  wave  velocity  may  be  a 
useful  index  of  liquefaction  potential  (Finn  1995a).  Shear  wave  velocity  is 
influenced  by  many  of  the  variables  that  influence  liquefaction,  such  as, 
soil  density,  confinement,  stress  history  and  geologic  age.  The  major 
advantage  of  using  shear  wave  velocity  as  an  index  of  liquefaction 
resistance  is  that  it  can  be  measured  in  soils  that  are  hard  to  sample,  such  as 
cohesionless  silt,  sand,  or  hard  to  penetrate  gravel. 

Currently,  seismic  methods  applied  to  the  problem  of  in  situ 
determination  of  liquefaction  potential  of  cohesionless  soils  are  limited  to 
the  determination  of  shear  wave  velocity.  Conventional  crosswell, 
downhole,  seismic  cone  penetrometer,  and  SASW  surface  wave  methods 
are  the  primary  seismic  tools  employed.  These  methods  are  one¬ 
dimensional,  yielding,  effectively,  velocity  values  for  a  vertical  sequence  of 
horizontal  ‘layers’  of  arbitrarily  chosen  thickness. 

Any  seismic  field  method  that  could  increase  the  accuracy  of  measured 
shear  wave  velocity,  or  expand  velocity  estimation  to  two  or  even  three 
dimensions  would  be  most  valuable.  Moreover,  seismic  methods  that 
could  produce  estimates  of  other  material  parameters,  such  as  density, 
porosity,  permeability,  frame  or  fluid  viscosity,  etc.,  would  be  extremely 
welcomed. 

Waveform  tomography  holds  the  promise  of  expansion  of  the  seismic 
liquefaction  index  beyond  shear  wave  velocity.  Currently,  waveform 
tomography  has  been  applied  to  the  simultaneous  evaluation  of  velocity 
and  density  for  an  acoustic  medium  (Mora  1987;  Sen  and  Stoffa  1991). 
While  problems  are  unresolved  relative  to  minimization  of  the  objective 
function  of  inversion  procedures,  the  potential  for  multi-parameter 
estimation  for  even  the  BISQ  model  is  promising. 

Field  tests  proposed  within  the  body  of  this  report  have  in  view 
advanced  data  processing  capability.  Two-dimensional,  crosswell, 
tomographic  data  gathering  exercises  are  proposed.  These  experiments 
include  borehole-to-surface  measurements  for  the  purpose  of  imaging 
shallow  targets.  In  addition,  two-dimensional,  surface  wave  experiments 
are  proposed.  One  exercise  will  be  conducted  for  collection  of  Rayleigh 
wave  data,  and  one  for  the  recording  of  Love  waves. 


While  awaiting  resolution  of  full  waveform  tomographic  inversion 
problems,  traveltime  attenuation  tomography  will  be  applied  to  measured 
crosswell  measurements  (including  borehole-to-surface  data). 

Tomographic  results  should  give  two-dimensional  distribution  estimates  of 
compressional  wave  velocity  and  attenuation,  soil  fabric  and  structure, 
shear  wave  velocity,  porosity,  and  possibly  permeability.  The  last  three 
items  will  be  extracted  from  compressional  wave  traveltime  and  amplitude 
data  by  means  of  empirical  relations  (e.g.,  Yamamoto  1994). 

It  is  experimentally  feasible  at  this  time  to  acquire  data  suitable  for  full 
waveform,  multi-parameter  tomography. 
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